;+
; NAME:
;       KDIST_SPECTRUM
;
; PURPOSE:
;       Calculates a probability density function for the distance to
;       a source given (l,b,v_spec).  Designed to be a modular insert
;       for prob_kidst and prob_grsmatch.
;
; CATEGORY:
;       distance-omnibus
;
; CALLING SEQUENCE:
;       pdf = KDIST_SPECTRUM(l, b, v_std, v_spec [,DVEC = dvec]
;                            [,ERRVLSR=errvlsr][,CONSTRAIN=constrain]
;                            [,/CLEM][,/IAU])
;
; INPUTS:
;       L         -- Galactic Longitude
;       B         -- Galactic Latitude
;       V_STD     -- Velocity values
;       V_SPEC    -- Velocity spectrum (used as weights)
;
; OPTIONAL INPUTS:
;       DVEC      -- Distance vector for calculation (Defualt: !MW
;                    parameters)
;       ERRVLSR   -- Estimated error in our determination of the VLSR
;                    (Default: value from !MW.ERRVLSR)
;
; KEYWORD PARAMETERS:
;       CLEM -- Use the Clemens (1985) rotation curve instead of the
;               Reid et al (2009) curve.
;       IAU  -- Use the IAU Standard rotation curve instead of the
;               Reid et al (2009) curve.
;
; OUTPUTS:
;       PDF       -- Distance pdf
;
; OPTIONAL OUTPUTS:
;       CONSTRAIN -- The constraint bit to be passed back up to the
;                    calling function.
;
; MODIFICATION HISTORY:
;       Created:  10/03/10, TPEB -- Initial Version to replace
;                                   identical code blocks in prob_kdist
;                                   and prob_grsmatch
;       Modified: 04/27/11, TPEB -- Clarified documentation, and
;                                   offloaded the conversion from
;                                   physical radial velocity to V_LSR
;                                   (as reported by telescopes) to a
;                                   new routine
;                                   (utils/vphys2vlsr.pro).
;       Modified: 04/04/12, TPEB -- Added pass-through keywords /CLEM
;                                   and /IAU for rotation curve
;                                   variants -- passed to VPHYS2VLSR.
;
;-

FUNCTION KDIST_SPECTRUM, l, b, v_std, v_spec, DVEC=dvec, ERRVLSR=errvlsr, $
                         CONSTRAIN=constrain, CLEM=clem, IAU=iau
  
  COMPILE_OPT OBSOLETE
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  ;; Set uniform probability in case of early return
  prob = dblarr(!mw.nbins)+1.d0/!mw.nbins
  constrain = 0b 
  
  ;; Check that v_spec is not null...
  IF TOTAL(v_spec) EQ 0 THEN RETURN, prob
  
  ;; If v_spec is nonzero, then calculate!
  IF n_elements(dvec) NE 0 THEN  d = dvec  ELSE BEGIN
     d = dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART
     dvec = d
  ENDELSE
  
  ;; Define !MW quantities into more useable variables
  IF n_elements(errvlsr) EQ 0 THEN errvlsr = !MW.ERRVLSR
                                ; Assume a 7 km/s error in the velocities
  
  
  ;; Use the routine VPHYS2VLSR to get the array of the VLSR that a
  ;; telescope would report (VLSR for each of the distances in the
  ;; dvec array).  This code block was moved to a new routine for
  ;; ease of use by other programs.
  vlsr = VPHYS2VLSR(l, b, d, CLEM=clem, IAU=iau)
  
  ;; Sum the probabilities, weighted by the v_spec.
  prob *= 0.d                              ;; Start w/ null prob 
  FOR q = 0L, n_elements(v_std)-1 DO BEGIN
     prob = prob + exp(-(vlsr-v_std[q])^2/(2*errvlsr^2))*v_spec[q]
  ENDFOR
  
  ;; Normalize & return
  prob = prob/total(prob)
  constrain = 1b
  RETURN, prob
  
END
